home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / MNBRAK.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  79 lines
  1. PROCEDURE mnbrak(VAR ax,bx,cx,fa,fb,fc: real);
  2. (* Programs using routine MNBRAK must supply an external
  3. function func(x:real):real for which a minimum is to be found *)
  4. LABEL 1;
  5. CONST
  6.    gold=1.618034;
  7.    glimit=100.0;
  8.    tiny=1.0e-20;
  9. VAR
  10.    ulim,u,r,q,fu,dum: real;
  11. FUNCTION max(a,b: real): real;
  12.    BEGIN
  13.       IF (a > b) THEN max := a ELSE max := b
  14.    END;
  15. FUNCTION sign(a,b: real): real;
  16.    BEGIN
  17.       IF (b > 0.0) THEN sign := abs(a) ELSE sign := -abs(a)
  18.    END;
  19. BEGIN
  20.    fa := func(ax);
  21.    fb := func(bx);
  22.    IF (fb > fa) THEN BEGIN
  23.       dum := ax;
  24.       ax := bx;
  25.       bx := dum;
  26.       dum := fb;
  27.       fb := fa;
  28.       fa := dum
  29.    END;
  30.    cx := bx+gold*(bx-ax);
  31.    fc := func(cx);
  32. 1:   IF (fb >= fc) THEN BEGIN
  33.       r := (bx-ax)*(fb-fc);
  34.       q := (bx-cx)*(fb-fa);
  35.       u := bx-((bx-cx)*q-(bx-ax)*r)/
  36.          (2.0*sign(max(abs(q-r),tiny),q-r));
  37.       ulim := bx+glimit*(cx-bx);
  38.       IF ((bx-u)*(u-cx) > 0.0) THEN BEGIN
  39.          fu := func(u);
  40.          IF (fu < fc) THEN BEGIN
  41.             ax := bx;
  42.             fa := fb;
  43.             bx := u;
  44.             fb := fu;
  45.             GOTO 1 END
  46.          ELSE IF (fu > fb) THEN BEGIN
  47.             cx := u;
  48.             fc := fu;
  49.             GOTO 1
  50.          END;
  51.          u := cx+gold*(cx-bx);
  52.          fu := func(u)
  53.       END ELSE IF  ((cx-u)*(u-ulim) > 0.0) THEN BEGIN
  54.          fu := func(u);
  55.          IF (fu < fc) THEN BEGIN
  56.             bx := cx;
  57.             cx := u;
  58.             u := cx+gold*(cx-bx);
  59.             fb := fc;
  60.             fc := fu;
  61.             fu := func(u)
  62.          END
  63.       END ELSE IF  ((u-ulim)*(ulim-cx) >= 0.0) THEN BEGIN
  64.          u := ulim;
  65.          fu := func(u)
  66.       END ELSE BEGIN
  67.          u := cx+gold*(cx-bx);
  68.          fu := func(u)
  69.       END;
  70.       ax := bx;
  71.       bx := cx;
  72.       cx := u;
  73.       fa := fb;
  74.       fb := fc;
  75.       fc := fu;
  76.       GOTO 1
  77.    END
  78. END;
  79.